Out-of-equilibrium bosons on a one-dimensional optical random lattice 
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We study the transport properties of a one-dimensional hard-core boson lattice gas coupled to 
two particle reservoirs at different chemical potentials generating a current flow through the system. 
In particular, the influence of random fluctuations of the underlying lattice on the stationary state 
properties is investigated. We show analytically that the steady-state density presents a linear 
proflle. The local steady-state current obeys the Fourier law j = — ft;(r)V/9 where r is a typical 
timescale of the lattice fluctuations and Vp the density gradient imposed by the reservoirs. 



Ultracold atoms/molecules trapped in optical lattices 
have become a very active field of both experimen- 
tal and theoretical research. For example, optical lat- 
tices can now be used to experimentally generate one- 
dimensional (ID) bosonic systems P, H Q which have 
been studied theoretically for many years 0, 0, 
In the large scattering-length limit and at low densi- 
ties, ultracold bosons effectively behave as impenetra- 
ble particles [3], namely as hard-core bosons, thus re- 
alising the Tonks-Girardeau model ^] . Experiments 
on such ID hard-core bosons have been performed with 
Rubidium atoms within both continuum and lat- 
tice contexts Q. In particular, theoretical and expri- 
mental considerations [ii, [13] show that a harmonically- 
trapped hard-core boson gas does not relax towards an 
equilibrium state even if the momentum distribution 
of an expanding-gas state ll| approaches that of non- 
interacting fermions Also worthy of mention is work 
on the influence of finite currents on the superfluid to 
Mott-insulator transition in a Bose-Hubbard model [3] ■ 
For a recent review of developments in ultracold gases 
and optical lattices, see In this letter we present the 
transport properties of a ID-lattice hard-core bosonic gas 
driven out of equilibrium by the interaction at its bound- 
aries with two external reservoirs that induce a particle- 
current flow through the system. We focus particularly 
on the influence of lattice fluctuations, which may be in- 
duced artificially or inherent to an experimental set-up, 
on the transport properties. 

The Hamiltonian of ID hard-core bosons is given by 
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where TV is the number of lattice sites, 6; , 6/ are the 
bosonic creation and annihilation operators at site I, 
ni = hfbi is the particle number operator, e the chemi- 
cal potential, and i/ the particle hopping rate across the 
Zth bond. The creation and annihilation operators sat- 
isfy the usual bosonic commutation relations on differ- 
ent sites, [bijlf,] — [bi,bi] — — 0, while the 
hard-core constraint is implemented by the additional 



conditions bf — (6^)^ = and {bi,b^} = 1 prevent- 
ing more than single occupancy of sites. Through the 
transformation b'^ = [a^ + ia'^)/2, where cr^^^ are the 
usual Pauli matrices, the hard-core boson Hamiltonian 
is exactly mapped onto the XX quantum spin chain 

which, in recent years, has been studied extensively in a 
nonequilibrium context [l5| . 

Our system is coupled to a quantum environment con- 
sisting of two independent reservoirs (at the left and the 
right of the system). The system- reservoir interaction is 
implemented via a repeated interaction scheme described 
as follows (l7|. At a given time a left (right) reservoir 



"particle" comes close to the lcft(right) boundary, inter- 
acts locally with the system during a short time inter- 
val T and goes away for ever, being replaced by another 
fresh reservoir particle. The process is repeated again 
and again. Let TLs be the system Hilbert space, Hn 
the Hilbert space of the nth environment copy with as- 
sociated Hamiltonian _ff„, and He ~ X^n* total 
environment Hamiltonian defined on f = 7i„ . We as- 
sume that initially, at i < 0, the system and the environ- 
ment are uncoupled and described by the factorized state 
p(0) — ps(0) ® Pe{0) where ps and pE are respectively 
the system and reservoir density matrices. By assump- 
tion, the reservoir state is factorized into pe{0) = 
where p„ is the initial state of the nth environment sub- 
system. At t > the interaction between the system 
and the reservoirs is switched on and the total Hamil- 
tonian on the time interval ](n — 1)t, nr] is given by 
= Hs + He + hY^^ where _ff j"^ is the interaction 
Hamiltonian between the nth copy of the environment 
and the system. The time evolution on that interval is 
Kn{T) = J/j"^ ® riNnfnl "^^^^ unitary evo- 



lution Ukir) = e"'^"" and U] 
The dynamics from < = to i nr is generated by the 
string operator U{nT) = Kn{T)Kn-i{T) . . . Ki{t) , and 
in the Schrodinger picture the total state evolves as 



p(nr) = ll{nT)p{0)ll\nT) ^ X„(r)p((n - 1)t)KI{i 
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The reduced density matrix associated with the sys- 
tem itself, ps{t) — TrE{p{t)}, evolves as psinr) = 

Tr„ l^uj"^[ps{{n - l)r) ® p„]C/j"^^| where Tr„ stands 

for the partial trace over the nth environment subsystem. 
Introducing the eigenstates of p„, {(/'^"^}j=i...dim{-H„} 
with eigenvalues Pj^\ such that Pn\<f>f^^) ~ pj"^ I'/'j"^), 

and the decomposition {\ip§) (8) |(?!>j"^)} of the joined nth 
system and subsystem Hilbert space 7is ® Tin, we ar- 
rive at the closed dynamics psinr) = K-n {ps{{n — 1)t)) 
for the system. Here /C„ is the completely positive map 
JCnX = YlTm^"^^ Pi^^^m^ym where the coefficients 
1/^, which are operators on Tis^ are defined through the 
decomposition of C/j"^ on {iV'f) ® \<f>Y^)}'- 
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Note that, for typographical convenience, we have sup- 
pressed here the dependence of the operators V on n. For 
identical subsystems, with p„ = pe Vn, we have = /C 
Vn. Iterating the previous recursive equation, we arrive 
at ps{nT) = /C" (P5(0)) . The Heisenberg equation of mo- 
tion is obtained through the adjoint action of K, defined 
through the trace-like scalar product (X, Y) = Tr{XY^} 
via {X^KY) = {IC^X,Y). Finally, for the time evolution 
of a system observable one has X{nT) = /Ct"X(0) . 

In our setup, the interaction with the reservoirs is 
given by the local left- and right-hopping Hamiltonians 

Hl = -tLibth + b+bL], Hr = -tR[b+bN-l + b^^.bRl 

describing the exchange of particles through the bound- 
aries between the reservoirs and the system. In the fol- 
lowing we set II = tR = Ie- The left and right reservoir 
particles are described by the one-particle density ma- 



trices pl,r. = (zL,fl) 



|0)|7o""(0| where 



zl,r are the left and right fugacities, |0), |1} the vacuum 
and one particle states and po: Pi their respective occu- 
pation probabilities, with pi+ pa = 1. 

The hard-core boson Hamiltonian ([T]) is canonically 
diagonalized by introducing a fermionic representation 
of the creation and annihilation operators &+, b [l6| : 



Ts] = e"^'f->^{bi+b+) 
Tsf ^te'^^"-^''^{bi~b+) 



(4) 
(5) 



where the F^s are Majorana real (Clifford) opera- 
tors satisfying F^''' = Fs and {F5",F5^} = 26ijSa0. 
In the Clifford representation, the Hamiltonian ^ 
takes the form Hs = {l/4)r s^TsTs where Ts^ = 
(rs},rs2, ■•■,rsjv)rsi, ...,Fs^) is a 27V-component op- 
erator and Ts ^ ^ 



written in terms of the tridiagonal matrix [Csjjfe = 
—i[e6jk + tj{Sjk+i + i5jfc-i)] with boundary conditions 
to ^ tN ~ 0- The time evolution of Fs, generated by Hs, 
is simply given by Tsit) = e-''^^TsiO) = Rs(i)rs(0), 
which defines the rotation matrix Rg. Its matrix ele- 
ments are expressed in terms of the spectral properties 



of Hs, see for the explicit forms. 



is a 2N X 2N matrix 



In the repeated interaction procedure, where the nth 
copy of the environment interacts with the system on 
time-interval ](n — 1)t, nr], the interaction Hamiltonian 
Hs + nj^^ + Hn is again of the form ((!]) with one 
more site, either on the left of the system if it is a left 
reservoir particle that interacts or on the right if it is a 
right reservoir particle. We introduce correspondingly a 
2(A^-|- l)-component Clifford vector operator F^"'' whose 
dynamics is governed by the new 2{N -I- 1) x 2(iV -I- 1) 
interaction matrix T^"> via rf">(T) = e"*^^'"'ri">. 
Exactly at time nr the expectation value of an 
observable Q, which has a given decomposition 
Q = FQ(r{">) on the F^">s, is given by (Q)(nT) = 

Trs,n |FQ(ri">)C/j">[p.((n - 1)t) ® p„]C/)"^^ 

Trs,„{FQ(r{"}(T))[p,((n-l)r)®p„]} . Due to the 
initial Gaussian factorized state describing the system 
and the environment, the reduced density matrix 
associated with the system itself remains a Gaussian 
state during the repeated interaction process. One has 
Ps{nT) = 2^^e-3rs^^^("^)^s .^viiere Ts{nT) is the 
time-evolved system interaction matrix. This implies 
that one may completely characterize the system by its 
two-point correlation matrix Gs = {—i^s'^s) + h where 
we have subtracted the trivial diagonal part by adding i 
to the definition. In the same way, the correlation matrix 
fjf"} = ^_jp{"}r{"}^ -(_ j completely characterizes the 
state of the system interacting with the nth environment 
copy. Using the somewhat shortened notation psn for 
ps{{n — l)r) ® pn, one has the dynamical equation 
G'{">(nT) = R'^-^^ {t)G^'^^ {{n - l)r)i?W(r)t with 



- Trs,n {-»ri">(r)F^;>(r)p5n} + i5k,k' (6) 
Gtia^ - 1)^) = Trs^n [-iT]^^t]^^ psn) + i5k,k' , (7) 



where we use the convention F^ for the fcth component 
of the vector F. The correlation matrix G^"''((n — 1)t) 
gives the two-point correlation just before the nth inter- 
action takes place, while the matrix G^"^(nT) contains 
the correlations after that interaction has taken place. 
Ordering the r{">^ = (r^,rt,) such that the first part, 
Te, is associated with the interacting part of the envi- 
ronment and the second part, F^, with the components 
of the system, one notices that the correlation matrix 
G^"j!, ((n — 1)t) takes a block-diagonal form [due to the 
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uncorrelated state — 1)t) ® 



G{">((n-l)r)= 



Gis 
Gs((n-l)r) 
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Decomposing the rotation matrix as R''^"^ 
Re Res 
RsE Rs 

namical equation for the system correlation matrix 



one arrives at the fundamental dy- 



Gsinr) = Rs{T)Gs{{n-l)T)RsHr)+RsE{T)GERsE\T) 

(9) 

Here Rs is a 2N x 2N square matrix and Re a 2 x 2 
one, while Res is a 2 x 2N rectangular matrix and Rse 
a 27V X 2 one. Notice that, for non-interacting dynamics, 
Res — 0, Rse = and the rotation matrix splits into 
a block-diagonal form where Rs — Rs and Re are the 
rotation matrices of the system and environment part re- 
spectively. Note also that in the dynamical equation ([9]), 
the bath properties enter only through the initial envi- 
ronment particle state, Ge, which stays constant in time 
since at each step of the repeated interaction procedure 
the bath particle is replaced by a fresh one. 

In the following, we concentrate mainly on the steady- 
state properties of the bosonic system, in particular its 
density profile and particle current. The hopping dynam- 
ics conserves particles and so one may define the cur- 
rent through the Heisenberg equation of motion hi = 
i[HN,ni] = J7j_i — Ji where Ji denotes the particle- 
current operator associated with the Ith bond and given 
by Ji = tiJk = iti [bibl^_^ - b^k+i] which is easily ex- 
pressed in terms of the Clifford operators F. 

If the lattice is free of any disorder then, since the parti- 
cles are identical, the hopping rates are uniform, i.e., ti — 
ts VZ. In this case, independent of the physical parame- 
ters ts/e, tE and r, the steady-state density profile is flat 
except at the sites directly in contact with the reservoir 
particles. The density in the flat region is the mean value 
set by the reservoirs n*i = p = {pL + Pr)/'2. V/ 7^ 1, TV [T1| 
while the boundary values are n\ = p — A{ts){pR — pl)/'^ 
and n*j^ = p + A{ts){pR — Pl)/'^ with a shift from the 
mean p depending on the strength ts and density dif- 
ference PR — PL- The steady-state current j* takes a 
constant value j* = —■^^{pr — Pl) independent of the 
system-size, which can be related to the fact that the to- 
tal current commutes with the system Hamiltonian. In 
this case there is no finite conductivity k and the system 
obviously does not obey Fourier's law. 

Next we consider the effect of fluctuations of the lattice 
parameters, which may well be induced artificially, on the 
steady-state properties. In particular, we take the limit 
of a strongly locahzed lattice gas, with ts/e,tE/£ 0, 
where the effect of the fluctuations is to enhance signifi- 
cantly the local hopping rates. The simplest case of this 
scenario is to consider that only one bond is activated 
during a given time r which is of the order of the inter- 
action time with the reservoirs. At each timestep r, a 



single bond is activated at random and locally the parti- 
cles are exchanged with a rate ts = 1/2, either within the 
system if the selected bond is a system one or with the 
reservoirs if the fluctuations act close to the boundaries. 
We introduce the expectation values of the density gradi- 
ent on the Ith bond: 61 = (71/+ 1) — and the associated 
current ji = (Ji) where ( . ) = Trs{ ■ Ps}- In the weak- 
hopping limit ts 0, the density gradient on bond I is 
changed only if the activated bond isl — l,lorl + l. From 
the dynamical equation ([9|), if the hopping is enhanced 
on bond I then Si is mapped to 6'i = cos t6i + sin rji and 
ji is mapped to jl — — sin tSi + cos r j; . Alternatively, if 
the activated bond is ^ ± 1, the updated density gradient 
satisfies 6[ = Si — ^(cost^j-i-i -I- sinrj;-|-i — Si±i) while 
the updated current is given by j; — cos(r/2)j; + Pi±i 
where the Pi±i oc F;F;-|-2 are proportional to correlations 
across two bonds. Under this dynamics, the system re- 
laxes exponentiall y to ward a current-carrying nonequilib- 
rium steady-state [18| . Nevertheless, we remark here that 
the periodicity of the dynamics implies a strong slowing- 
down of the relaxation to the steady-state in the neigh- 
borhood of r = n2n. For even values of n the dynamical 
generator maps to the identity and for odd values it maps 
to a reflection dynamics which loses its relaxation prop- 
erties; below we avoid these special r values. 

In the steady state, the P terms vanish on average and 
the set of dynamical equations for the gradient density 
and current closes. At any given time, the last update on 
bond I has probability 1/3 to have resulted from an up- 
date on bond I and probability 1/3(1/3) to have resulted 
from an update on bond I — l{l + 1). Consequently, the 
steady-state averaged (denoted by a star) gradient obeys 



(cost — l)Si + sinrj;* = rj^r) 



(10) 



where ry(r) is a constant independent of the bond in- 
dex I. Since the steady-state current — j* is con- 
stant in space it also follows from pU]) that the gra- 
dient density is site-independent and we can thus omit 
the /-subscripts. The steady-state current satisfies j* = 
■1 [— sinT(5* + cos rj*] + | cos ^j* which is equivalent to 
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3 — cos T — 2 cos ^ 



6* 



Mr)S* (11) 



and defines the conductivity coefficient k{t). We now 
determine the constant 77(t) by considering the bound- 
ary conditions. Remembering that the densities on the 
reservoir sites are fixed and that the boundary currents 
jo = JN = due to the repeated interaction scheme, one 
sees that an update on bond (between left reservoir and 
boundary site) gives (5g = ^(1 + cosr)(5o whereas an up- 
date on bond 1 gives (5q = <5o — ^(cosr^i -I- sinrji — 5i). 
The steady-state average must therefore obey 



1 



cos T 



cosT(5* + sinrj* - 5* 







(12) 
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FIG. 1: Normalized steady-state density gradient as a func- 
tion of the interaction time with a delta time-distribution on 
the left and an exponential one with mean To on the right. 
The full lines correspond to the analytical curves while the 
crosses are obtained numerically with a time average of the 
density gradient. 



giving 77(t) = (cost — 1)6q . By symmetry, at the right- 
boundary we find Sn = 6q. Noting that the reservoir 
density diflFerence Ap = pR — pi^ = 25q + [N — 1)(5*, one 
finally gets for the bulk steady-state density gradient 



5* 



Ap 



iV+ 1 + 7(r) 



(13) 



with the finite-size shift function (extrapolation length) 
7(t) — 2 ^^^J^^ k{t). This analytical expression is com- 
pared with numerical simulation data in Fig. 1 and the 
agreement is seen to be excellent. 

So far, we have considered the somewhat unphysical 
situation where the enhancement of a local hopping rate 
always stays precisely for a time r. This hypothesis leads 
to the trigonometric form of the conductivity k{t) and 
shift function 7(t). A more reasonable assumption would 
be to draw the duration of a local fluctuation from a prob- 
ability distribution /(r). In that case, one may average 
the dynamical equations for the current and the gradient 
density over the time-distribution /(r) which basically 
leads to replacing the trigonometric functions cost, sinT 
and cos ^ by their expectations under /. For example, 
for an exponential distribution of interaction timescales 
with mean Tq, one gets the shift fimction 7 = (2/to)k 

with conductivity k = ^^^l^^'^2) '^hich is again in agree- 
ment with the numerical results, see Fig. 1. In the limit 
of short-time interaction Tq — > 0, conductivity diverges 
as T~^ while the density gradient goes as t^ leading to a 
linear vanishing of the steady-state current j* ~ Tq. 

At finite, but small, hopping rate ti = ts along the 
optical lattice, we observe numerically that the linear 
density profile survives. However, the gradient density 
is strongly attenuated by a function, of the bulk hop- 
ping rate ts, whose asymptotic behaviour is {Nts)~^ for 
large lattice sizes. Consequently, for large systems the 



steady-state gradient density behaves as 5* ~ 1/7V^ in- 
stead of the former 1/N behaviour. At the same time, 
the steady-state current j* ^ leading to a Hnear 

divergence of the conductivity coefficient k with system 
size. This implies that the classical transport properties 
of the hard-core boson gas are induced by the optical 
lattice fluctuations only for sufficiently small Nts values. 

In conclusion, we have derived analytical expressions 
for the steady-state conductivity and density profile of 
a nonequilibrium hard-core boson model. For a perfect 
optical lattice, due to the integrability of the model, the 
transport properties are anomalous with an infinite con- 
ductivity coefficient. On the other hand, when fluctu- 
ations of the underlying lattice are present and when 
Nts ^ 1, the classical Fourier law is recovered, with a a 
linear density profile and a finite conductivity coefficient. 
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